// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H

#include "./Tridiagonalization.h"

namespace Eigen {

/** \eigenvalues_module \ingroup Eigenvalues_Module
 *
 *
 * \class GeneralizedSelfAdjointEigenSolver
 *
 * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
 *
 * \tparam _MatrixType the type of the matrix of which we are computing the
 * eigendecomposition; this is expected to be an instantiation of the Matrix
 * class template.
 *
 * This class solves the generalized eigenvalue problem
 * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
 * selfadjoint and the matrix \f$ B \f$ should be positive definite.
 *
 * Only the \b lower \b triangular \b part of the input matrix is referenced.
 *
 * Call the function compute() to compute the eigenvalues and eigenvectors of
 * a given matrix. Alternatively, you can use the
 * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
 * constructor which computes the eigenvalues and eigenvectors at construction time.
 * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
 * and eigenvectors() functions.
 *
 * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
 * contains an example of the typical use of this class.
 *
 * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
 */
template<typename _MatrixType>
class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
{
	typedef SelfAdjointEigenSolver<_MatrixType> Base;

  public:
	typedef _MatrixType MatrixType;

	/** \brief Default constructor for fixed-size matrices.
	 *
	 * The default constructor is useful in cases in which the user intends to
	 * perform decompositions via compute(). This constructor
	 * can only be used if \p _MatrixType is a fixed-size matrix; use
	 * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
	 */
	GeneralizedSelfAdjointEigenSolver()
		: Base()
	{
	}

	/** \brief Constructor, pre-allocates memory for dynamic-size matrices.
	 *
	 * \param [in]  size  Positive integer, size of the matrix whose
	 * eigenvalues and eigenvectors will be computed.
	 *
	 * This constructor is useful for dynamic-size matrices, when the user
	 * intends to perform decompositions via compute(). The \p size
	 * parameter is only used as a hint. It is not an error to give a wrong
	 * \p size, but it may impair performance.
	 *
	 * \sa compute() for an example
	 */
	explicit GeneralizedSelfAdjointEigenSolver(Index size)
		: Base(size)
	{
	}

	/** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
	 *
	 * \param[in]  matA  Selfadjoint matrix in matrix pencil.
	 *                   Only the lower triangular part of the matrix is referenced.
	 * \param[in]  matB  Positive-definite matrix in matrix pencil.
	 *                   Only the lower triangular part of the matrix is referenced.
	 * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
	 *                     Default is #ComputeEigenvectors|#Ax_lBx.
	 *
	 * This constructor calls compute(const MatrixType&, const MatrixType&, int)
	 * to compute the eigenvalues and (if requested) the eigenvectors of the
	 * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
	 * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
	 * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
	 * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
	 * \a options contains ComputeEigenvectors.
	 *
	 * In addition, the two following variants can be solved via \p options:
	 * - \c ABx_lx: \f$ ABx = \lambda x \f$
	 * - \c BAx_lx: \f$ BAx = \lambda x \f$
	 *
	 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
	 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
	 *
	 * \sa compute(const MatrixType&, const MatrixType&, int)
	 */
	GeneralizedSelfAdjointEigenSolver(const MatrixType& matA,
									  const MatrixType& matB,
									  int options = ComputeEigenvectors | Ax_lBx)
		: Base(matA.cols())
	{
		compute(matA, matB, options);
	}

	/** \brief Computes generalized eigendecomposition of given matrix pencil.
	 *
	 * \param[in]  matA  Selfadjoint matrix in matrix pencil.
	 *                   Only the lower triangular part of the matrix is referenced.
	 * \param[in]  matB  Positive-definite matrix in matrix pencil.
	 *                   Only the lower triangular part of the matrix is referenced.
	 * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
	 *                     Default is #ComputeEigenvectors|#Ax_lBx.
	 *
	 * \returns    Reference to \c *this
	 *
	 * According to \p options, this function computes eigenvalues and (if requested)
	 * the eigenvectors of one of the following three generalized eigenproblems:
	 * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
	 * - \c ABx_lx: \f$ ABx = \lambda x \f$
	 * - \c BAx_lx: \f$ BAx = \lambda x \f$
	 * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
	 * matrix \f$ B \f$.
	 * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
	 *
	 * The eigenvalues() function can be used to retrieve
	 * the eigenvalues. If \p options contains ComputeEigenvectors, then the
	 * eigenvectors are also computed and can be retrieved by calling
	 * eigenvectors().
	 *
	 * The implementation uses LLT to compute the Cholesky decomposition
	 * \f$ B = LL^* \f$ and computes the classical eigendecomposition
	 * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
	 * and of \f$ L^{*} A L \f$ otherwise. This solves the
	 * generalized eigenproblem, because any solution of the generalized
	 * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
	 * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
	 * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
	 * can be made for the two other variants.
	 *
	 * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
	 * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
	 *
	 * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
	 */
	GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA,
											   const MatrixType& matB,
											   int options = ComputeEigenvectors | Ax_lBx);

  protected:
};

template<typename MatrixType>
GeneralizedSelfAdjointEigenSolver<MatrixType>&
GeneralizedSelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matA, const MatrixType& matB, int options)
{
	eigen_assert(matA.cols() == matA.rows() && matB.rows() == matA.rows() && matB.cols() == matB.rows());
	eigen_assert((options & ~(EigVecMask | GenEigMask)) == 0 && (options & EigVecMask) != EigVecMask &&
				 ((options & GenEigMask) == 0 || (options & GenEigMask) == Ax_lBx || (options & GenEigMask) == ABx_lx ||
				  (options & GenEigMask) == BAx_lx) &&
				 "invalid option parameter");

	bool computeEigVecs = ((options & EigVecMask) == 0) || ((options & EigVecMask) == ComputeEigenvectors);

	// Compute the cholesky decomposition of matB = L L' = U'U
	LLT<MatrixType> cholB(matB);

	int type = (options & GenEigMask);
	if (type == 0)
		type = Ax_lBx;

	if (type == Ax_lBx) {
		// compute C = inv(L) A inv(L')
		MatrixType matC = matA.template selfadjointView<Lower>();
		cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
		cholB.matrixU().template solveInPlace<OnTheRight>(matC);

		Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);

		// transform back the eigen vectors: evecs = inv(U) * evecs
		if (computeEigVecs)
			cholB.matrixU().solveInPlace(Base::m_eivec);
	} else if (type == ABx_lx) {
		// compute C = L' A L
		MatrixType matC = matA.template selfadjointView<Lower>();
		matC = matC * cholB.matrixL();
		matC = cholB.matrixU() * matC;

		Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);

		// transform back the eigen vectors: evecs = inv(U) * evecs
		if (computeEigVecs)
			cholB.matrixU().solveInPlace(Base::m_eivec);
	} else if (type == BAx_lx) {
		// compute C = L' A L
		MatrixType matC = matA.template selfadjointView<Lower>();
		matC = matC * cholB.matrixL();
		matC = cholB.matrixU() * matC;

		Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);

		// transform back the eigen vectors: evecs = L * evecs
		if (computeEigVecs)
			Base::m_eivec = cholB.matrixL() * Base::m_eivec;
	}

	return *this;
}

} // end namespace Eigen

#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
